Indicators
The following indicators are used in this report
- Surface and bottom temperature anomalies (FVCOM, OISST)
- Surface and bottom salinity anomalies (FVCOM)
- Maine Coastal Current Index (FVCOM)
- Species-based lobster predator index (NOAA Trawl Survey)
- Continuous Plankton Recordings
# Indicators
FVCOM <- read_csv(here::here("Indicators/FVCOM_stat_area_temps.csv")) %>% filter(Year >= 1980)
OISST <- read_csv(here::here("Indicators/OISST_stat_area_anoms.csv")) %>% filter(Date <= as.Date("2020-12-31"))
sal <- read_csv(here::here("Indicators/FVCOM_stat_area_sal.csv")) %>% filter(Year >= 1980)
maineCC <- read_csv(here::here("Indicators/mcc_pca_pc1_2.csv")) %>%
rename("Year" = yr, "Month" = mon)
species_index <- read_rds(here::here("Indicators/nmfs_trawl_lobster_predator_indices.rdata")) %>%
rename("Year" = year)
mcc_turnoff_subset <- read_csv("/Users/mdzaugis/Documents/EcosystemIndicators/Indicators/mcc_turnoff_subset.csv")
cprData <- read_csv(here::here("Indicators/cpr_focal_pca_timeseries_period_1961-20017.csv")) %>%
mutate(`First Mode` = `First Mode`*-1) %>%
rename("Year" = year,
"FirstMode" = `First Mode`,
"SecondMode" = `Second Mode`) %>%
select(Year, FirstMode, SecondMode)
Temperature
Source: FVCOM NECOFS Monthly Means Thredds Link Methods:
- Surface and bottom layer (1 m above bathymetry)
- Regridded to regular 0.1 deg grid
- Averaged over NOAA statistical area
- Baseline climatology 1990-2020
yrOISST <- OISST %>%
mutate(Year = lubridate::year(Date)) %>%
group_by(Year, stat_area) %>%
summarise(temp_var = var(temp),
temp = mean(temp),
.groups = "drop") %>%
mutate(stat_area = as.factor(stat_area))
yrFVCOM <- FVCOM %>%
group_by(Year, stat_area) %>%
summarise(sur_temp = mean(sur_temp_anom),
bot_temp = mean(bot_temp_anom),
sur_temp_var = var(sur_temp_anom),
bot_temp_var = var(bot_temp_anom),
.groups="drop") %>%
mutate(stat_area = as.factor(stat_area))
tPlot1 <- yrOISST %>%
ggplot() +
geom_line(aes(Year, temp, col = stat_area)) +
theme_bw() +
labs(y = "SST anomalies (deg C)")
tPlot2 <- yrFVCOM %>%
ggplot() +
geom_line(aes(Year, bot_temp, col = stat_area)) +
theme_bw() +
labs(y = "Bottom T anomalies (deg C)")
tPlot1/tPlot2

Salinity
Source: FVCOM NECOFS Monthly Means Thredds Link Methods:
- Surface and bottom layer (1 m above bathymetry)
- Regridded to regular 0.1 deg grid
- Averaged over NOAA statistical area for each year
- Baseline climatology 1990-2017 (last year of data)
yrSal <- sal %>%
group_by(Year, stat_area) %>%
summarise(sur_sal = mean(sur_sal_anom),
bot_sal = mean(bot_sal_anom),
sur_sal_var = var(sur_sal_anom),
bot_sal_var = var(bot_sal_anom), .groups = "drop") %>%
mutate(stat_area = as.factor(stat_area))
sPlot1 <- yrSal %>%
ggplot() +
geom_line(aes(Year, sur_sal, col = stat_area)) +
theme_bw() +
labs(y = "SSS anomalies (deg C)")
sPlot2 <- yrSal %>%
ggplot() +
geom_line(aes(Year, bot_sal, col = stat_area)) +
theme_bw() +
labs(y = "Bottom S anomalies (deg C)")
sPlot1/sPlot2

Maine Coastal Current Index
Source: FVCOM NECOFS Monthly Means Thredds Link Methods:
- Surface layer
- Crop to Maine Coastal Current interest area
- Regridded to regular 0.1 deg grid
- Averaged over NOAA statistical area for each year
- See MCC_index_report.Rmd for details
summerMCC <- maineCC %>%
filter(Month %in% c(6,7,8,9)) %>%
group_by(Year) %>%
summarise(MCCPC1 = mean(PC1))
yrMCC <- maineCC %>%
group_by(Year) %>%
summarise(MCCPC1 = mean(PC1))
MCCPC_plot <- yrMCC %>%
mutate(c = ifelse(MCCPC1 > 0, 1, 2)) %>%
ggplot() + geom_col(aes(x = Year, y = MCCPC1, fill = as.factor(c))) + theme_bw() + theme(legend.position = "none") + labs(x = "Year", y = "Maine Coastal Current PC1")
MCCPC_maps <- mcc_turnoff_subset %>%
mutate(Year = lubridate::year(as.Date(Date)),
Month = lubridate::month(as.Date(Date))) %>%
filter(Year %in% c(1980, 1987, 2000, 2011, 2012),
Month %in% c(6,7,8,9)) %>%
group_by(Year, lat, lon) %>%
summarise(u = mean(u),
v = mean(v), .groups = "drop") %>%
mutate(vel = sqrt(u^2+v^2),
PC = if_else(Year == 2010, "negative PC1", if_else(Year ==2013, "positive PC1", "neutral PC1"))) %>%
ggplot() + geom_sf(data= ne_us, fill = "grey") +
geom_segment(aes(x = lon, y = lat, xend=lon+u, yend=lat+v, color = vel),
arrow = arrow(angle = 10, length = unit(0.05, "inches"), type = "closed")) +
scale_color_viridis_c() + theme_bw() +
coord_sf(datum = "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0") +
facet_wrap(~Year, nrow = 1)+
theme(panel.background = element_blank(), panel.grid = element_blank(),
axis.title = element_blank(), axis.text = element_blank(),
axis.ticks = element_blank())
MCCPC_plot / MCCPC_maps

Species Based Predator Index
Source: NOAA NE Fisheries Trawl Survey Methods:
- Filtered to 15 lobster predators (ASMFC 2020)
- Cropped to NOAA stat areas 511, 512, 513
- Stratum abundance: abundance per trawl summed over stat area
- Stratified abundance: abundance / km2 multiplied by the area of the strata
yrSpecies_index <-
species_index %>%
filter(season == "Both",
`stratum id` != 'Strata 511-513') %>%
dplyr::select(Year, `stratum id`, "predators" = `stratified biomass (kg)`)
yrSpecies_index %>%
ggplot() +
geom_line(aes(Year, predators, color = `stratum id` )) + theme_bw()

CPR data
Source: Continuous Plankton Recording
- First Mode = smaller zooplankton species explains ~ 50% of variance
- Second Mode = Calanus explains ~ 26% of variance
cprData %>%
pivot_longer(cols = c(FirstMode, SecondMode),
names_to = "Mode", values_to = "values") %>%
ggplot() +
geom_line(aes(Year, values, col = Mode)) +
theme_bw() +
scale_color_discrete(label = c("Small spp", "Calanus"))

# area weighted
# regional time series
# area weighted average
Indicators PCA
allIndex <- left_join(yrFVCOM, yrSal, by = c("Year", "stat_area")) %>%
left_join(yrOISST, by = c("Year", "stat_area")) %>%
left_join(yrMCC, by = c("Year" = "Year")) %>%
left_join(yrSpecies_index, by = c("Year" = "Year",
"stat_area" = "stratum id")) %>%
left_join(cprData, by = "Year") %>%
na.omit() %>%
select(-sur_temp, -sur_sal, -sur_sal_var, -bot_sal_var, -sur_temp_var, -bot_temp_var, -temp_var, -temp)
Previous analysis shows surface and bottom temperature and surface and bottom salinity load similarly in a PCA. To reduce variables only bottom temperature and salinity are used in the following PCA.
- Data: Bottom temp, bottom salinity, MCC, Predator Index, CPR data mode one and two
- Method: Principal components analysis
- Years: 1990-2016
PCA Summary
For each stat area, the first two principal components have eigenvalues greater than one and were retained for further analysis. PC1 and PC2 explain 34.8% and 24.0% of the variability, respectively, in stat area 511, 30.5% and 23.2% in stat area 512, and 33.2% and 23.3% in stat area 513.
# PCA by stat area
indicator_pca <- allIndex %>%
nest(data = c(-stat_area),
Year = c(Year)) %>%
mutate(data = map(data, ~column_to_rownames(.x, var = "Year")),
pca = map(data, ~prcomp(.x, scale. = TRUE, center = TRUE)),
pca_index = map(pca, ~tibble("PC1" = .x$x[,1],
"PC2" = .x$x[,2],
"PC3" = .x$x[,3])))
index_pca <- indicator_pca %>%
unnest(c(pca_index, Year, data)) %>%
mutate(PC2 = if_else(stat_area == "511", PC2, PC2),
PC1 = if_else(stat_area == "512", PC1, PC1),
PC3 = if_else(stat_area == "513", PC3, PC3))
PCA511 <- data.frame(summary(indicator_pca$pca[[1]])$importance) %>%
rownames_to_column("Result") %>% mutate(stat_area = "511")
PCA512 <- data.frame(summary(indicator_pca$pca[[2]])$importance) %>%
rownames_to_column("Result") %>% mutate(stat_area = "512")
PCA513 <- data.frame(summary(indicator_pca$pca[[3]])$importance) %>%
rownames_to_column("Result") %>% mutate(stat_area = "513")
PCA_table <- bind_rows(PCA511,PCA512,PCA513)
knitr::kable(PCA_table)
| Standard deviation |
1.444888 |
1.200781 |
0.9993327 |
0.7902938 |
0.7644566 |
0.5126402 |
511 |
| Proportion of Variance |
0.347950 |
0.240310 |
0.1664400 |
0.1040900 |
0.0974000 |
0.0438000 |
511 |
| Cumulative Proportion |
0.347950 |
0.588260 |
0.7547100 |
0.8588000 |
0.9562000 |
1.0000000 |
511 |
| Standard deviation |
1.353820 |
1.178900 |
1.0157298 |
0.8498775 |
0.8141508 |
0.6004372 |
512 |
| Proportion of Variance |
0.305470 |
0.231630 |
0.1719500 |
0.1203800 |
0.1104700 |
0.0600900 |
512 |
| Cumulative Proportion |
0.305470 |
0.537110 |
0.7090600 |
0.8294400 |
0.9399100 |
1.0000000 |
512 |
| Standard deviation |
1.410562 |
1.182894 |
1.0209041 |
0.8574840 |
0.7567640 |
0.5107456 |
513 |
| Proportion of Variance |
0.331610 |
0.233210 |
0.1737100 |
0.1225500 |
0.0954500 |
0.0434800 |
513 |
| Cumulative Proportion |
0.331610 |
0.564820 |
0.7385300 |
0.8610700 |
0.9565200 |
1.0000000 |
513 |
scree1 <- fviz_eig(indicator_pca$pca[[1]], ylab = "% variance explained") + scale_y_continuous(limits = c(0,60))
scree2 <- fviz_eig(indicator_pca$pca[[2]], ylab = "% variance explained") + scale_y_continuous(limits = c(0,60))
scree3 <- fviz_eig(indicator_pca$pca[[3]], ylab = "% variance explained") + scale_y_continuous(limits = c(0,60))
#scree1 / scree2 / scree3
PC time series
Below are the time series for PC1, PC2, and PC3. PC1 and PC2 time series are highly correlated between stat areas.
PC1plot <- index_pca %>%
ggplot() +
geom_line(aes(Year, PC1, color = stat_area))
PC2plot <- index_pca %>%
ggplot() +
geom_line(aes(Year, PC2, color = stat_area))
PC3plot <- index_pca %>%
ggplot() +
geom_line(aes(Year, PC3, color = stat_area))
pc1cors <- index_pca %>%
select(stat_area, PC1, PC2, PC3, Year) %>%
pivot_wider(names_from = stat_area, values_from = c(PC1, PC2, PC3)) %>%
summarise("511_512" = corrr::correlate(PC1_511, PC1_512)[[2]],
"511_513" = corrr::correlate(PC1_511, PC1_513)[[2]],
"512_513" = corrr::correlate(PC1_512, PC1_513)[[2]])
pc2cors <- index_pca %>%
select(stat_area, PC1, PC2, PC3, Year) %>%
pivot_wider(names_from = stat_area, values_from = c(PC1, PC2, PC3)) %>%
summarise("511_512" = corrr::correlate(PC2_511, PC2_512)[[2]],
"511_513" = corrr::correlate(PC2_511, PC2_513)[[2]],
"512_513" = corrr::correlate(PC2_512, PC2_513)[[2]])
pc3cors <- index_pca %>%
select(stat_area, PC1, PC2, PC3, Year) %>%
pivot_wider(names_from = stat_area, values_from = c(PC1, PC2, PC3)) %>%
summarise("511_512" = corrr::correlate(PC3_511, PC3_512)[[2]],
"511_513" = corrr::correlate(PC3_511, PC1_513)[[2]],
"512_513" = corrr::correlate(PC3_512, PC3_513)[[2]])
PC1plot / PC2plot / PC3plot

Loadings plot
Indicators load similarly for each stat area. The Maine coastal current, first mode (small zooplankton), and second mode (Calanus) tend to load together. The lobster predator index has the most variable relationship to PC1 and PC2 with respect to stat area. In stat area 512 and 513, predators loads with bottom temperature. In area 511, predators load separated from the rest of the indicators and has a stronger influence on PC2 than PC1. Bottom salinity and the plankton FirstMode (Calanus) tend to load similarly for each stat area.
Interpretation:
- PC1: Bottom salinity and Maine Coastal Current load directly opposite and are the largest contributors to PC1
- PC1: Contribution of bottom temperature increases from stat area 511 to 513
- PC2: Small zooplankton and Calanus group together and load opposite to bottom temperature and predators
| bot_temp |
0.393 |
0.458 |
0.537 |
-0.475 |
-0.418 |
-0.348 |
0.214 |
0.055 |
0.220 |
| bot_sal |
0.610 |
0.612 |
0.546 |
0.013 |
0.103 |
0.232 |
0.199 |
0.038 |
-0.144 |
| MCCPC1 |
-0.575 |
-0.583 |
-0.462 |
-0.083 |
-0.107 |
-0.260 |
-0.139 |
0.190 |
0.486 |
| predators |
-0.264 |
0.132 |
0.395 |
-0.298 |
-0.306 |
-0.447 |
0.762 |
0.813 |
0.408 |
| FirstMode |
0.223 |
0.235 |
0.209 |
0.611 |
0.599 |
0.581 |
0.008 |
-0.056 |
0.146 |
| SecondMode |
-0.153 |
-0.058 |
-0.009 |
0.552 |
0.593 |
0.469 |
0.560 |
0.543 |
0.712 |
Chronolgical cluster
The cluster plot groups years that are most similar based on PC1 and PC2. For this analysis stat areas were grouped together.
clusfun <- function(x){
wss <- (nrow(x)-1)*sum(apply(x,2,var)) # get sum of squares
for (i in 2:12) wss[i] <- sum(kmeans(x,
centers=i)$withinss)
return(wss)
}
allIndex_ca <- allIndex %>%
nest(data = -stat_area) %>%
mutate(data = map(data, ~column_to_rownames(.x, var = "Year")),
data = map(data, ~scale(.x)),
wss = map(data, ~clusfun(.x)))
# look for break in plot like a scree plot
kmeans_scree <- allIndex_ca %>%
unnest(wss) %>%
add_column("x" = rep(c(1:12),3)) %>%
ggplot() +
geom_line(aes(x = x, y = wss, col = stat_area)) +
labs(x = "Number of Clusters",
y ="Within groups sum of squares")
# K-Means Cluster Analysis
fit <- allIndex_ca %>%
mutate(fit = map(data, ~kmeans(.x, 3))) # 4 cluster solution
(fviz_cluster(object = fit$fit[[1]], data = fit$data[[1]]) +
theme_bw() + biplot511) /
(fviz_cluster(object = fit$fit[[2]], data = fit$data[[2]]) +
theme_bw() + biplot512) /
(fviz_cluster(object = fit$fit[[3]], data = fit$data[[3]]) +
theme_bw() + biplot513)

Breakpoint Analysis
The breakpoint analysis estimates a change in the linear relationship in the data. The location of the break indicates there may be a difference in the relationship of the variable before and after that point. We find that breakpoints change depending on the variable.
Breakpoint of indicators
# breapoint function
bp_analysis <- function(x){
mod <- lm(values ~ Year, data = x)
o <- tryCatch(segmented::segmented(mod, seg.Z = ~Year, psi = c(2005)), # need to estimate bp
error = function(cond){cond})
}
set.seed(25)
# Breakpoint by stat area
indicator_breakpoint <- allIndex %>%
pivot_longer(cols = c(bot_temp, bot_sal, MCCPC1, predators, FirstMode, SecondMode),
names_to = "Indicator", values_to = "values") %>%
nest(data = c(-stat_area, -Indicator)) %>%
mutate(seg = purrr::map(data, ~bp_analysis(.x)),
bp = purrr::map(seg, ~data.frame(.x$psi)))
# get fitted values and plot over data
bp_plots <- list()
for(i in 1:length(indicator_breakpoint$data)){
name_bp <- paste0(indicator_breakpoint$stat_area[[i]], "_", indicator_breakpoint$Indicator[[i]])
# get the fitted data
fitted_df <- fitted(indicator_breakpoint$seg[[i]])
if(is.null(fitted_df)){
"check error log"
} else {
seg_model <- data.frame(Year = indicator_breakpoint$data[[i]]$Year, values = fitted_df)
# plot the fitted model
bp_plots[[name_bp]] <- indicator_breakpoint$data[[i]] %>%
ggplot() + geom_line(aes(Year, values)) +
geom_line(data = seg_model, aes(x = Year, y = values), col = "dark red")
}
}
# extract the breakpoints
break_years <- indicator_breakpoint %>%
unnest(bp) %>%
select(stat_area, Indicator, Est.) %>%
mutate(Indicator = factor(Indicator, levels = c("MCCPC1", "predators", "bot_sal", "bot_temp", "FirstMode", "SecondMode")))
# plot break points
ggplot(break_years) +
geom_point(aes(Est., Indicator, color = stat_area, shape = stat_area), size = 3) +
scale_color_grey(name = "Stat area") +
scale_shape(name = "Stat area") +
labs(x = "Estimated breakpoint")

Breakpoint of PC1 and PC2
set.seed(8)
# Breakpoint by stat area
pc1_breakpoint <- index_pca %>%
select(stat_area, Year, PC1, PC2) %>%
pivot_longer(cols = c(PC1, PC2),
names_to = "Indicator", values_to = "values") %>%
nest(data = c(-stat_area, -Indicator))%>%
mutate(seg = purrr::map(data, ~bp_analysis(.x)),
bp = purrr::map(seg, ~data.frame(.x$psi))) %>%
unnest(bp) %>%
select(stat_area, Indicator, Est., St.Err)
knitr::kable(pc1_breakpoint)
| 511 |
PC1 |
2004.646 |
1.557871 |
| 511 |
PC2 |
1994.000 |
5.878710 |
| 512 |
PC1 |
2004.729 |
1.780807 |
| 512 |
PC2 |
2009.000 |
9.367909 |
| 513 |
PC1 |
1997.193 |
2.662064 |
| 513 |
PC2 |
1998.000 |
4.618772 |
Lobster Data
ALSI <- read_csv(here::here("Biological_data/SettlementIndex.csv"))
ME_landings <- read_csv(here::here("Biological_data/ME_lob_landings_1950_2019.csv")) %>%
add_column("stat_area" = "State Wide")
ASFMCindicies <- read_csv(here::here("Biological_data/GOMGBK_indices.csv")) %>%
rename("lob_index" = Index, "name" = Survey)
GOMindex <- ASFMCindicies %>%
filter(name %in% c("MeFQ2", "MeMQ2")) %>%
add_column("stat_area" = "State Wide")
NEFSCindex <- ASFMCindicies %>%
filter(name %in% c("NefscFQ2", "NefscMQ2")) %>%
add_column("stat_area" = "State Wide")
sublegal <- readxl::read_excel(here::here("Biological_data/RawData.xlsx")) %>%
rename("Year" = year, "Month" = month, "stat_area" = `stat area`)
MEDMR_trawl <- read_csv(paste0(gmRi::box_path("Res_Data", "Maine_NH_Trawl"), "MaineDMR_Trawl_Survey_Catch_Data_2021-05-14.csv")) %>%
filter(Common_Name == "Lobster American")
NEFSC_meIndex <- read_csv(here::here("Biological_data/nmfs_trawl_lobster_indices.csv")) %>%
rename("Year" = est_year, "stat_area" = lobster_strata) %>%
mutate(stat_area = as.factor(stat_area))
ALSI index
Source: American Lobster Settlement Index data portal Sites:
- Jonesport, Length: 2002-2018, stat area: 511
- Mt. Desert Island, Length: 2000-2018, stat area: 512
- Outer Penobscot Bay, Length: 2000-2018, stat area: 512
- Mid-coast, Length: 1989-2018, stat area: 513
- Casco Bay, Length: 2000-2018, stat area: 513
- York, Length: 2000-2018, stat area: 513
Methods:
- Sites grouped by NOAA stat area
- Time series cropped to shortest length
- Averaged by stat area
statKey <- data.frame("Location" = c('Jonesport',
'Mt. Desert Island',
'Outer Pen Bay',
'Mid-coast',
'Casco Bay',
'York'), "stat_area" = c(511,512,512,513,513,513))
ALSI <- ALSI %>%
left_join(statKey, by = "Location") %>%
filter(!is.na(stat_area),
Year >= 2000) %>%
group_by(stat_area, Year) %>%
summarise(lob_index = mean(Yoy_density), .groups = "drop") %>%
mutate(stat_area = as.factor(stat_area),
name = "ALSI") %>%
na.omit()
ALSI_cors <- ALSI %>%
left_join(index_pca, by = c("Year", "stat_area")) %>%
group_by(stat_area, name) %>%
summarise(corPC1 = corrr::correlate(lob_index, PC1)[[2]],
corPC2 = corrr::correlate(lob_index, PC2)[[2]],
corPC3 = corrr::correlate(lob_index, PC3)[[2]])
ALSI_plot <- ALSI %>%
ggplot() +
geom_line(aes(Year, lob_index, col = as.factor(stat_area))) +
scale_color_discrete(name = "Stat area") +
labs(y = "ALSI yoy density")
bp_analysis <- function(x, psi){
mod <- lm(lob_index ~ Year, data = x)
o <- tryCatch(segmented::segmented(mod, seg.Z = ~Year, psi = psi), # need to estimate bp
error = function(cond){cond})
}
ALSI_bp <- ALSI %>%
nest(data = -c(stat_area, name)) %>%
mutate(mod = map(data, ~bp_analysis(.x, c(2005))),
bp1 = purrr::map(mod, ~.x$psi[[2]]))
ALSI_plot

Sublegal index
Source: Ventless trap survey
Calculate stratified means
- Calculate catch per unit effort for each site for each year
- Multiply cpue by the depth strata area factor
- Group by stat area, sum the outputs
strata_area <- tibble("stat_area" = c(511,512,513),
"1" = c(122, 566, 315),
"2" = c(82, 395, 338),
"3" = c(92, 420, 198)) %>%
pivot_longer(cols = c("1", "2", "3"), names_to = "depth stratum", values_to = "strata_area") %>%
group_by(stat_area) %>%
mutate(`depth stratum` = as.numeric(`depth stratum`),
total = sum(strata_area))
# average number of lobsters per trap haul at each site
u <- sublegal %>%
group_by(Year, `site ID`, `effort ID`, stat_area, `depth stratum`) %>%
summarise(n_lob = sum(`sample ID` != 0), .groups = "drop") %>%
group_by(Year, `site ID`, stat_area, `depth stratum`) %>%
summarise(n_lob_u = mean(n_lob), .groups = "drop")
# average number of lobsters per trap haul at each depth stratum within a stat area
v <- sublegal %>%
group_by(stat_area, `depth stratum`, Year, `site ID`, `effort ID`) %>%
summarise(n_lob = sum(`sample ID` != 0), .groups = "drop") %>%
group_by(stat_area, `depth stratum`, Year) %>%
summarise(n_lob_v = mean(n_lob), .groups = "drop")
# choose the relevant v which corresponds to the depth strata that u is in ("stat_area", "depth stratum", "Year"). Each site will have one w value
w <- left_join(v, u, by = c("stat_area", "depth stratum", "Year")) %>%
mutate(w = (n_lob_v-n_lob_u)^2)
#sum of all w within the same depth strata in a stat area
x <- w %>%
group_by(`depth stratum`, stat_area, Year) %>%
summarise(x = sum(w), .groups = "drop")
#number of sites per depth stratum within a given stat area
y <- sublegal %>%
group_by(stat_area, `depth stratum`, Year) %>%
summarise(site_ID = unique(`site ID`)) %>%
summarise(y = n(), .groups = "drop")
#if done correctly, each stat area per year should have three z values, one for each depth stratum
z <- left_join(x, y, by = c("stat_area", "depth stratum", "Year")) %>%
mutate(z = x/(y-1))
# Calculate stat area variance
sublegal_vari <- left_join(z, strata_area, by = c("stat_area", "depth stratum")) %>%
mutate(a = 1/(total^2),
b = strata_area*(strata_area-y)*(z/y)) %>%
group_by(stat_area, Year, a) %>%
summarise(stat_sum = sum(b), .groups = "drop") %>%
mutate(vari = a*stat_sum,
sd = sqrt(vari))
strata_area_factor <- tibble("stat_area" = c(511,512,513),
"1" = c(0.412162162, 0.409847936, 0.370152761),
"2" = c(0.277027027, 0.28602462, 0.397179788),
"3" = c(0.310810811, 0.304127444, 0.23266745)) %>%
pivot_longer(cols = c("1", "2", "3"), names_to = "depth stratum", values_to = "strata_scale") %>%
mutate(`depth stratum` = as.numeric(`depth stratum`))
sublegal_cpue <- sublegal %>%
group_by(`trip ID`, `trip date`, Year, Month, `site ID`,
`depth stratum`, stat_area,
`soak nights`, depth, `depth unit`,
`latitude (dd)`, `longitude (dd)`, `effort ID`) %>%
mutate(lob_count = sum(`sample ID` != 0),
lob_effort = lob_count/`soak nights`) %>%
group_by(`trip ID`, `trip date`, Year, Month,
`site ID`, `depth stratum`, stat_area,
`soak nights`, depth, `depth unit`,
`latitude (dd)`, `longitude (dd)`) %>%
summarise(cpue = sum(lob_effort)/sum(!is.na(unique(`effort ID`))), .groups = "drop") %>%
left_join(strata_area_factor, by = c("stat_area", "depth stratum")) %>%
mutate(stratified_cpue = cpue*strata_scale,
stat_area = as.factor(stat_area)) %>%
group_by(stat_area, Year) %>%
summarise(lob_index = sum(stratified_cpue), .groups = "drop") %>%
mutate(name = "sublegal_cpue")
sub_leagal_cors <- sublegal_cpue %>%
left_join(index_pca, by = c("Year", "stat_area")) %>%
group_by(stat_area, name) %>%
summarise(corPC1 = corrr::correlate(lob_index, PC1)[[2]],
corPC2 = corrr::correlate(lob_index, PC2)[[2]],
corPC3 = corrr::correlate(lob_index, PC3)[[2]])
sublegal_plot <- sublegal_cpue %>%
ggplot() +
geom_line(aes(Year, lob_index, col = as.factor(stat_area))) +
scale_color_discrete(name = "Stat area") +
labs(y = "Sublegal lobster cpue")
sublegal_var_plot <- sublegal_vari %>%
ggplot() +
geom_line(aes(Year, vari, col = as.factor(stat_area))) +
scale_color_discrete(name = "Stat area") +
labs(y = "Sublegal lobster variance")
sublegal_bp <- sublegal_cpue %>%
nest(data = -c(stat_area, name)) %>%
mutate(mod = map(data, ~bp_analysis(.x, c(2016))),
bp1 = purrr::map(mod, ~.x$psi[[2]]))
sublegal_plot

ME-NH Trawl Survey
Source: ME-NH Trawl Survey
Calculate stratified means
- Calculate catch per unit effort for each site for each year
- Multiply cpue by the depth strata area factor
- Group by stat area, sum the outputs
stat_area_trawl_key <- tibble("stat_area" = c(511, 512, 512, 513, 513),
"Region" = c(1,2,3,4,5))
DMR_strata_area <- tibble("Stratum" = c("1", "2", "3", "4"),
"1" = c(253.27, 214.22, 227.35, 225.65),
"2" = c(279.63, 191.23, 211.66, 263.49),
"3" = c(259.62, 262.90, 280.03, 183.69),
"4" = c(205.30, 206.12, 310.49, 170.72),
"5" = c(138.54, 220.49, 365.04, 196.11)) %>%
pivot_longer(cols = c("1", "2", "3", "4", "5"), names_to = "Region", values_to = "strata_area") %>%
group_by(Region) %>%
mutate(Stratum = as.numeric(Stratum),
total = sum(strata_area),
Region = as.numeric(Region)) %>%
left_join(stat_area_trawl_key) %>%
group_by(stat_area, Stratum) %>%
summarise(strata_area = sum(strata_area),
total = sum(total))
# average number of lobsters per tow
u <- MEDMR_trawl %>%
left_join(stat_area_trawl_key) %>%
group_by(Season, Year, Tow_Number, stat_area, Stratum) %>%
summarise(n_lob_u = sum(Expanded_Weight_kg, na.rm = TRUE), .groups = "drop")
# average number of lobsters per trap haul at each depth stratum within a stat area
v <- MEDMR_trawl %>%
left_join(stat_area_trawl_key) %>%
group_by(stat_area, Stratum, Year, Season) %>%
summarise(n_lob_v = mean(Expanded_Weight_kg, na.rm = TRUE), .groups = "drop")
# choose the relevant v which corresponds to the depth strata that u is in ("stat_area", "depth stratum", "Year"). Each site will have one w value
w <- left_join(v, u, by = c("stat_area", "Stratum", "Year", "Season")) %>%
mutate(w = (n_lob_v-n_lob_u)^2)
#sum of all w within the same depth strata in a stat area
x <- w %>%
group_by(Stratum, stat_area, Year, Season) %>%
summarise(x = sum(w), .groups = "drop")
#number of sites per depth stratum within a given stat area
y <- MEDMR_trawl %>%
left_join(stat_area_trawl_key) %>%
group_by(stat_area, Stratum, Year, Season) %>%
summarise(Tow_number = unique(Tow_Number, na.rm = TRUE)) %>%
summarise(y = n(), .groups = "drop")
#if done correctly, each stat area per year should have three z values, one for each depth stratum
z <- left_join(x, y, by = c("stat_area", "Stratum", "Year", "Season")) %>%
mutate(z = x/(y-1))
# Calculate stat area variance
MEDMR_vari <- left_join(z, DMR_strata_area, by = c("stat_area", "Stratum")) %>%
mutate(a = 1/(total^2),
b = strata_area*(strata_area-y)*(z/y)) %>%
group_by(stat_area, Year, a, Season) %>%
summarise(stat_sum = sum(b), .groups = "drop") %>%
mutate(vari = a*stat_sum,
sd = sqrt(vari))
medmr_vari_plot <- MEDMR_vari %>%
ggplot() +
geom_line(aes(Year, vari, col = as.factor(stat_area))) +
facet_wrap(~Season)
MEDMR_cpue <- MEDMR_trawl %>%
left_join(stat_area_trawl_key) %>%
group_by(Year, Season, stat_area, Stratum)%>%
mutate(tows=n_distinct(Tow_Number))%>%
group_by(Year, Season, tows, stat_area, Stratum) %>%
summarise(weight = sum(Expanded_Weight_kg,na.rm=TRUE),
catch=sum(Expanded_Catch,na.rm=TRUE))%>%
mutate(weight_tow = weight/tows,
catch_tow = catch/tows) %>%
left_join(DMR_strata_area) %>%
mutate(stratified_wpue = weight_tow*(strata_area/total),
stratified_cpue = catch_tow*(strata_area/total)) %>%
group_by(stat_area, Year, Season) %>%
summarise(lob_cpue = sum(stratified_cpue),
lob_index = sum(stratified_wpue),
.groups = "drop") %>%
mutate(stat_area = as.factor(stat_area),
name = "ME-NH_trawl")
MEDMR_cors <- MEDMR_cpue %>%
left_join(index_pca, by = c("Year", "stat_area")) %>%
group_by(stat_area, Season) %>%
summarise(corPC1 = corrr::correlate(lob_index, PC1)[[2]],
corPC2 = corrr::correlate(lob_index, PC2)[[2]],
corPC3 = corrr::correlate(lob_index, PC3)[[2]])
MEDMR_cpue %>%
ggplot() +
geom_line(aes(Year, lob_index, col = as.factor(stat_area))) +
scale_color_discrete(name = "Stat area") +
labs(y = "ME-NH Trawl lobster cpue") +
facet_wrap(~Season)

MEDMR_bp <- MEDMR_cpue %>%
nest(data = -c(stat_area, Season, name)) %>%
mutate(mod = map(data, ~bp_analysis(.x, c(2015))),
bp1 = purrr::map(mod, ~.x$psi[[2]]))
ME yearly landings
Yearly Maine lobster landings in pounds from 1950-2020.
ME_landings_cors <- ME_landings %>%
left_join(index_pca, by = c("Year", "stat_area")) %>%
group_by(stat_area) %>%
summarise(corPC1 = corrr::correlate(Pounds, PC1)[[2]],
corPC2 = corrr::correlate(Pounds, PC2)[[2]],
corPC3 = corrr::correlate(Pounds, PC3)[[2]])
ME_pounds <- ME_landings %>%
mutate(name = "ME_landings") %>%
select(Year, "lob_index" = Pounds, name, stat_area)
ME_pounds %>%
ggplot() +
geom_line(aes(Year, lob_index)) +
labs(y = "Maine lobster landings (pounds)")

ME_pounds_bp <- ME_pounds %>%
nest(data = -c(stat_area, name)) %>%
mutate(mod = map(data, ~bp_analysis(.x, c(1990, 2015))),
bp1 = purrr::map(mod, ~.x$psi[[2]]))
Nefsc index
ASFMC lobster abundance index based on the NE fisheries trawl survey. Time spans 1978-2018 and indices are split into seasons and sex.
NEFSC_cors <- NEFSCindex %>%
left_join(index_pca, by = c("Year", "stat_area")) %>%
group_by(stat_area, name) %>%
summarise(corPC1 = corrr::correlate(lob_index, PC1)[[2]],
corPC2 = corrr::correlate(lob_index, PC2)[[2]],
corPC3 = corrr::correlate(lob_index, PC3)[[2]], .groups = "drop")
NEFSCindex %>%
ggplot() +
geom_line(aes(Year, lob_index)) +
facet_wrap(~name) +
labs(y = "Nefsc trawl lob abundance index")

NEFSCindex_bp <- NEFSCindex %>%
nest(data = -c(stat_area, name)) %>%
mutate(mod = map(data, ~bp_analysis(.x, c(1990, 2015))),
bp1 = purrr::map(mod, ~.x$psi[[2]]))
NMFS stat area 511-513
Lobster biomass and abundance for statistical areas 511-513
NEFSC_meIndex <- NEFSC_meIndex %>% rename("lob_index" = `stratified biomass (kg)`)
NEFSC_meIndex_cors <- NEFSC_meIndex %>%
left_join(index_pca, by = c("Year", "stat_area")) %>%
group_by(stat_area) %>%
summarise(corPC1 = corrr::correlate(lob_index, PC1)[[2]],
corPC2 = corrr::correlate(lob_index, PC2)[[2]],
corPC3 = corrr::correlate(lob_index, PC3)[[2]], .groups = "drop")
nmfsBIO <- NEFSC_meIndex %>%
ggplot() +
geom_line(aes(Year, lob_index, col = as.factor(stat_area))) +
facet_wrap(~season) +
labs(y = "Lobster biomass") +
scale_color_discrete(name = "Stat area")
nmfsBIO_var <- NEFSC_meIndex %>%
ggplot() +
geom_line(aes(Year, bio_var, col = as.factor(stat_area))) +
facet_wrap(~season) +
labs(y = "variance") +
scale_color_discrete(name = "Stat area")
nmfsNUM <- NEFSC_meIndex %>%
ggplot() +
geom_line(aes(Year, `stratified abundance`, col = as.factor(stat_area))) +
facet_wrap(~season) +
labs(y = "Lobster abundance") +
scale_color_discrete(name = "Stat area")
nmfsNUM_var <- NEFSC_meIndex %>%
ggplot() +
geom_line(aes(Year, abun_var, col = as.factor(stat_area))) +
facet_wrap(~season) +
labs(y = "variance") +
scale_color_discrete(name = "Stat area")
nmfsBIO / nmfsNUM

NEFSC_meIndex_bp <- NEFSC_meIndex %>%
nest(data = -c(stat_area, season)) %>%
mutate(mod = map(data, ~bp_analysis(.x, c(1990, 2015))),
bp1 = purrr::map(mod, ~.x$psi[[1,2]]),
bp2 = purrr::map(mod, ~.x$psi[[2,2]]))
Biological data analysis
xyPlot <- function(x){
if(is.character(x$stat_area)){
x <- select(x, -stat_area)
}
x %>%
left_join(index_pca) %>%
na.omit() %>%
ggplot() +
geom_point(aes(PC1, lob_index, color = stat_area)) +
geom_smooth(aes(PC1, lob_index, color = stat_area), method = "gam", se = FALSE) +
facet_wrap(~name)
}
stepwiseAIC <- function(x){
if(is.character(x$stat_area)){
x <- select(x, -stat_area)
}
x %>%
left_join(index_pca) %>%
na.omit() %>%
select(stat_area, lob_index, PC1, PC2, PC3, name) %>%
group_by(stat_area, name) %>%
nest() %>%
mutate(aic = purrr::map(data, ~MASS::stepAIC(object = lm(lob_index ~ PC1 + PC2 + PC3, data = .x), direction = "both", trace = 0)),
stp = purrr::map(aic, broom::glance),
index = purrr::map(aic, ~as.character(.x$call$formula)),
model = paste(index[[1]][[2]], index[[1]][[1]], index[[1]][[3]])) %>%
select(stat_area, stp, model, name) %>%
unnest(c(stp, model))
}
aicModel <- function(x){
if(is.character(x$stat_area)){
x <- select(x, -stat_area)
}
x %>%
left_join(index_pca) %>%
na.omit() %>%
select(stat_area, lob_index, PC1, PC2, PC3, name) %>%
group_by(stat_area, name) %>%
nest() %>%
mutate(aic = purrr::map(data, ~MASS::stepAIC(object = lm(lob_index ~ PC1 + PC2 + PC3, data = .x), direction = "both", trace = 0)),
stp = purrr::map(aic, broom::glance),
index = purrr::map(aic, ~as.character(.x$call$formula)),
model = paste(index[[1]][[2]], index[[1]][[1]], index[[1]][[3]])) %>%
select(stat_area, stp, model, name) %>%
unnest(c(stp, model))
}
aicModelGAM <- function(x, k){
if(is.character(x$stat_area)){
x <- select(x, -stat_area)
}
x %>%
left_join(index_pca) %>%
na.omit() %>%
select(stat_area, lob_index, PC1, PC2, PC3, name) %>%
group_by(stat_area, name) %>%
nest() %>%
mutate(gams = purrr::map(data, ~mgcv::gam(lob_index ~ s(PC1, k=k) + s(PC2, k=k) + s(PC3, k=k),
data = .x, select = TRUE, method="REML")),
tidied = purrr::map(gams, broom::tidy)) %>%
select(stat_area, name, tidied) %>%
unnest(c(tidied))
}
aicModel_indicators <- function(x){
if(is.character(x$stat_area)){
x <- select(x, -stat_area)
}
x %>%
left_join(index_pca) %>%
na.omit() %>%
select(stat_area, lob_index, bot_temp, bot_sal, MCCPC1, predators, FirstMode, SecondMode, name) %>%
group_by(stat_area, name) %>%
nest() %>%
mutate(aic = purrr::map(data, ~MASS::stepAIC(object = lm(lob_index ~ bot_temp + bot_sal + MCCPC1 + predators + FirstMode + SecondMode, data = .x), direction = "both", trace = 0)),
stp = purrr::map(aic, broom::glance),
index = purrr::map(aic, ~as.character(.x$call$formula)),
model = paste(index[[1]][[2]], index[[1]][[1]], index[[1]][[3]])) %>%
select(stat_area, stp, model, name) %>%
unnest(c(stp, model))
}
piecewise_reg <- function(x, bps){
x <- x %>% unnest(c(data, bp1)) %>%
select(stat_area, name, Year, lob_index, bp1)
if(is.character(x$stat_area)){
x <- select(x, -stat_area)
}
x <- x %>%
left_join(index_pca) %>%
na.omit() %>%
select(stat_area, lob_index, PC1, PC2, PC3, name, bp1, Year) %>%
mutate(kd = if_else(Year>bp1, 1, 0)) %>%
group_by(stat_area, name) %>%
nest() %>%
mutate(aic = purrr::map(data, ~MASS::stepAIC(object = lm(lob_index ~ PC1*(DK) +
PC2*(DK) +
PC3*(DK), data = .x),
direction = "both", trace = 0)),
stp = purrr::map(aic, broom::glance),
index = purrr::map(aic, ~as.character(.x$call$formula)),
model = paste(index[[1]][[2]], index[[1]][[1]], index[[1]][[3]]))
x <- x %>%
unnest(data)
}
##################################################
ALSI index
The ALSI index is not significantly correlated with PC1 or PC
Correlation table
ALSI_cors %>%
na.omit() %>%
knitr::kable()
| 511 |
ALSI |
-0.4194353 |
0.3072101 |
0.1653596 |
| 512 |
ALSI |
-0.2646401 |
0.3738632 |
0.1615255 |
| 513 |
ALSI |
-0.4993881 |
0.6493318 |
0.5116346 |
Scatter plot
xyPlot(ALSI)

AIC lm model selection table
aicSel <- aicModel(ALSI)
aicSelGAM <- aicModelGAM(ALSI, 5)
aicSel2 <- aicModel_indicators(ALSI)
knitr::kable(aicSel)
| 511 |
0.1759260 |
0.1125357 |
0.3832675 |
2.775282 |
0.1196292 |
1 |
-5.825488 |
17.65098 |
19.77513 |
1.909621 |
13 |
15 |
lob_index ~ PC1 |
ALSI |
| 512 |
0.1397737 |
0.0824252 |
0.3313509 |
2.437271 |
0.1393292 |
1 |
-4.280251 |
14.56050 |
17.06014 |
1.646901 |
15 |
17 |
lob_index ~ PC2 |
ALSI |
| 513 |
0.6492657 |
0.5991608 |
0.3174624 |
12.958130 |
0.0006529 |
2 |
-2.965900 |
13.93180 |
17.26465 |
1.410953 |
14 |
17 |
lob_index ~ PC2 + PC3 |
ALSI |
knitr::kable(aicSelGAM)
| 511 |
ALSI |
s(PC1) |
1.4882941 |
4 |
1.4041080 |
0.0498928 |
| 511 |
ALSI |
s(PC2) |
0.0000041 |
4 |
0.0000002 |
0.7939962 |
| 511 |
ALSI |
s(PC3) |
0.0877091 |
4 |
0.0240327 |
0.3139331 |
| 512 |
ALSI |
s(PC1) |
0.0000062 |
4 |
0.0000010 |
0.4875469 |
| 512 |
ALSI |
s(PC2) |
0.5897069 |
4 |
0.3593165 |
0.1385407 |
| 512 |
ALSI |
s(PC3) |
0.0000041 |
4 |
0.0000005 |
0.5923790 |
| 513 |
ALSI |
s(PC1) |
0.0000028 |
4 |
0.0000004 |
0.5527103 |
| 513 |
ALSI |
s(PC2) |
0.9359963 |
4 |
3.6557366 |
0.0012759 |
| 513 |
ALSI |
s(PC3) |
0.8910265 |
4 |
2.0441191 |
0.0086693 |
Subleagal index
Correlation table
sub_leagal_cors %>%
na.omit() %>%
knitr::kable()
| 511 |
sublegal_cpue |
0.4932694 |
-0.5993328 |
0.5726005 |
| 512 |
sublegal_cpue |
0.5766013 |
-0.6320931 |
0.0502201 |
| 513 |
sublegal_cpue |
0.4611801 |
-0.5554480 |
-0.3723648 |
Scatter plot
xyPlot(sublegal_cpue)

AIC lm model selection table
aicSelGAM <- aicModelGAM(sublegal_cpue, 4)
aicSel <- aicModel(sublegal_cpue)
aicSel2 <- aicModel_indicators(sublegal_cpue)
knitr::kable(aicSel)
| 511 |
0.5910509 |
0.4888137 |
984.1015 |
5.781169 |
0.0279690 |
2 |
-89.66585 |
187.3317 |
188.9233 |
7747646 |
8 |
11 |
lob_index ~ PC2 + PC3 |
sublegal_cpue |
| 512 |
0.3995417 |
0.3328241 |
3207.2477 |
5.988552 |
0.0369289 |
1 |
-103.30949 |
212.6190 |
213.8127 |
92577938 |
9 |
11 |
lob_index ~ PC2 |
sublegal_cpue |
| 513 |
0.3085225 |
0.2316917 |
1464.6883 |
4.015608 |
0.0760731 |
1 |
-94.68801 |
195.3760 |
196.5697 |
19307808 |
9 |
11 |
lob_index ~ PC2 |
sublegal_cpue |
knitr::kable(aicSelGAM)
| 511 |
sublegal_cpue |
s(PC1) |
0.0000503 |
3 |
0.0000070 |
0.6157636 |
| 511 |
sublegal_cpue |
s(PC2) |
0.8202627 |
3 |
1.5162891 |
0.0440829 |
| 511 |
sublegal_cpue |
s(PC3) |
0.7952929 |
3 |
1.2948368 |
0.0557396 |
| 512 |
sublegal_cpue |
s(PC1) |
0.8877120 |
3 |
0.5732628 |
0.1778300 |
| 512 |
sublegal_cpue |
s(PC2) |
0.6330516 |
3 |
0.5747109 |
0.1049775 |
| 512 |
sublegal_cpue |
s(PC3) |
0.0000518 |
3 |
0.0000023 |
0.9114541 |
| 513 |
sublegal_cpue |
s(PC1) |
0.0001354 |
3 |
0.0000317 |
0.4395255 |
| 513 |
sublegal_cpue |
s(PC2) |
0.8145738 |
3 |
1.4640482 |
0.0414786 |
| 513 |
sublegal_cpue |
s(PC3) |
0.8303740 |
3 |
0.4817565 |
0.2221681 |
ME yearly landings
Correlation table
ME_landings_cors %>%
na.omit() %>%
knitr::kable()
Scatter plot
xyPlot(ME_pounds)

AIC lm model selection table
aicSelGAM <- aicModelGAM(ME_pounds, 5)
aicSel <- aicModel(ME_pounds)
aicSel2 <- aicModel_indicators(ME_pounds)
knitr::kable(aicSel)
| 511 |
0.2484089 |
0.2028579 |
32420331 |
5.453427 |
0.0089887 |
2 |
-672.1102 |
1352.220 |
1358.555 |
3.468557e+16 |
33 |
36 |
lob_index ~ PC1 + PC2 |
ME_landings |
| 512 |
0.1596119 |
0.1348946 |
33774123 |
6.457500 |
0.0157796 |
1 |
-674.1203 |
1354.241 |
1358.991 |
3.878351e+16 |
34 |
36 |
lob_index ~ PC2 |
ME_landings |
| 513 |
0.3718964 |
0.3338295 |
29637548 |
9.769553 |
0.0004651 |
2 |
-668.8795 |
1345.759 |
1352.093 |
2.898668e+16 |
33 |
36 |
lob_index ~ PC1 + PC2 |
ME_landings |
knitr::kable(aicSelGAM)
| 511 |
ME_landings |
s(PC1) |
1.1419314 |
4 |
0.6751062 |
0.1046332 |
| 511 |
ME_landings |
s(PC2) |
1.5806096 |
4 |
2.1875071 |
0.0071152 |
| 511 |
ME_landings |
s(PC3) |
0.0008035 |
4 |
0.0001255 |
0.4743764 |
| 512 |
ME_landings |
s(PC1) |
0.9048502 |
4 |
0.3531980 |
0.2155047 |
| 512 |
ME_landings |
s(PC2) |
0.8566724 |
4 |
1.4797510 |
0.0125743 |
| 512 |
ME_landings |
s(PC3) |
0.0003620 |
4 |
0.0000200 |
0.7621572 |
| 513 |
ME_landings |
s(PC1) |
0.8235815 |
4 |
1.1663465 |
0.0192908 |
| 513 |
ME_landings |
s(PC2) |
1.7420407 |
4 |
4.1600693 |
0.0003216 |
| 513 |
ME_landings |
s(PC3) |
0.6312029 |
4 |
0.2467791 |
0.2130456 |
ME Trawl index
Correlation table
MEDMR_cors %>%
na.omit() %>%
knitr::kable()
| 511 |
Fall |
0.6682315 |
-0.4816462 |
0.0686523 |
| 511 |
Spring |
0.6491385 |
-0.6148421 |
0.0777621 |
| 512 |
Fall |
0.2990573 |
-0.4904639 |
-0.3677374 |
| 512 |
Spring |
0.8260627 |
-0.7212692 |
-0.1826128 |
| 513 |
Fall |
0.2835740 |
-0.1419138 |
-0.6338938 |
| 513 |
Spring |
0.3547642 |
-0.2634394 |
-0.6122392 |
Scatter plot
xyPlot(MEDMR_cpue)

AIC lm model selection table
aicSel <- aicModel(MEDMR_cpue)
aicSelGAM <- aicModelGAM(MEDMR_cpue, 5)
aicSel2 <- aicModel_indicators(MEDMR_cpue)
knitr::kable(aicSel)
| 511 |
0.2520418 |
0.2279141 |
11.30559 |
10.446168 |
0.0029114 |
1 |
-125.8282 |
257.6564 |
262.1459 |
3962.305 |
31 |
33 |
lob_index ~ PC1 |
ME-NH_trawl |
| 512 |
0.3246583 |
0.2796355 |
15.32223 |
7.210977 |
0.0027724 |
2 |
-135.3194 |
278.6388 |
284.6249 |
7043.123 |
30 |
33 |
lob_index ~ PC2 + PC3 |
ME-NH_trawl |
| 513 |
0.3834734 |
0.3635855 |
22.56166 |
19.281692 |
0.0001219 |
1 |
-148.6297 |
303.2594 |
307.7489 |
15779.887 |
31 |
33 |
lob_index ~ PC3 |
ME-NH_trawl |
knitr::kable(aicSelGAM)
| 511 |
ME-NH_trawl |
s(PC1) |
1.5546299 |
4 |
1.9841589 |
0.0078687 |
| 511 |
ME-NH_trawl |
s(PC2) |
0.3831062 |
4 |
0.1549088 |
0.1922815 |
| 511 |
ME-NH_trawl |
s(PC3) |
0.0000867 |
4 |
0.0000094 |
0.5528045 |
| 512 |
ME-NH_trawl |
s(PC1) |
0.2418640 |
4 |
0.0703663 |
0.3080260 |
| 512 |
ME-NH_trawl |
s(PC2) |
0.9026573 |
4 |
2.3138608 |
0.0022853 |
| 512 |
ME-NH_trawl |
s(PC3) |
0.5777374 |
4 |
0.3419736 |
0.1333209 |
| 513 |
ME-NH_trawl |
s(PC1) |
0.8969429 |
4 |
0.3981576 |
0.1538907 |
| 513 |
ME-NH_trawl |
s(PC2) |
1.3321819 |
4 |
0.9360308 |
0.0568851 |
| 513 |
ME-NH_trawl |
s(PC3) |
0.9534783 |
4 |
5.1131060 |
0.0000568 |
NEFSC trawl index
Correlation table
NEFSC_cors %>%
filter(name %in% c("NefscFQ2", "NefscFQ4", "NefscMQ2", "NefscMQ4")) %>%
na.omit() %>%
knitr::kable()
Scatter plot
xyPlot(NEFSCindex)

AIC lm model selection table
aicSelGAM <- aicModelGAM(NEFSCindex, 5)
aicSel <- aicModel(NEFSCindex)
aicSel2 <- aicModel_indicators(NEFSCindex)
knitr::kable(aicSel)
| 511 |
0.3144294 |
0.2728796 |
1.966978 |
7.567543 |
0.0019718 |
2 |
-73.86951 |
155.7390 |
162.0731 |
127.67701 |
33 |
36 |
lob_index ~ PC1 + PC2 |
NefscFQ2 |
| 512 |
0.2476817 |
0.2020866 |
2.060507 |
5.432205 |
0.0091333 |
2 |
-75.54186 |
159.0837 |
165.4178 |
140.10775 |
33 |
36 |
lob_index ~ PC1 + PC2 |
NefscFQ2 |
| 513 |
0.3387701 |
0.2986955 |
1.931744 |
8.453499 |
0.0010860 |
2 |
-73.21882 |
154.4376 |
160.7717 |
123.14393 |
33 |
36 |
lob_index ~ PC1 + PC2 |
NefscFQ2 |
| 511 |
0.3780106 |
0.3403142 |
1.255375 |
10.027782 |
0.0003957 |
2 |
-57.70321 |
123.4064 |
129.7405 |
52.00687 |
33 |
36 |
lob_index ~ PC1 + PC2 |
NefscMQ2 |
| 512 |
0.3079663 |
0.2660249 |
1.324175 |
7.342770 |
0.0023020 |
2 |
-59.62402 |
127.2480 |
133.5821 |
57.86354 |
33 |
36 |
lob_index ~ PC1 + PC2 |
NefscMQ2 |
| 513 |
0.4109079 |
0.3752053 |
1.221725 |
11.509201 |
0.0001614 |
2 |
-56.72508 |
121.4502 |
127.7842 |
49.25621 |
33 |
36 |
lob_index ~ PC1 + PC2 |
NefscMQ2 |
knitr::kable(aicSelGAM)
| 511 |
NefscFQ2 |
s(PC1) |
1.2876032 |
4 |
1.4506571 |
0.0191078 |
| 511 |
NefscFQ2 |
s(PC2) |
1.5880929 |
4 |
2.2550748 |
0.0063438 |
| 511 |
NefscFQ2 |
s(PC3) |
0.0001103 |
4 |
0.0000148 |
0.5213610 |
| 512 |
NefscFQ2 |
s(PC1) |
0.7669306 |
4 |
0.8215864 |
0.0460963 |
| 512 |
NefscFQ2 |
s(PC2) |
0.8480125 |
4 |
1.3946331 |
0.0148886 |
| 512 |
NefscFQ2 |
s(PC3) |
0.0000609 |
4 |
0.0000080 |
0.5078377 |
| 513 |
NefscFQ2 |
s(PC1) |
0.8374938 |
4 |
1.2882553 |
0.0182335 |
| 513 |
NefscFQ2 |
s(PC2) |
0.9070660 |
4 |
2.4384850 |
0.0023881 |
| 513 |
NefscFQ2 |
s(PC3) |
0.0000246 |
4 |
0.0000039 |
0.4821666 |
| 511 |
NefscMQ2 |
s(PC1) |
1.4701779 |
4 |
2.0387306 |
0.0078542 |
| 511 |
NefscMQ2 |
s(PC2) |
0.9197827 |
4 |
2.8571796 |
0.0011273 |
| 511 |
NefscMQ2 |
s(PC3) |
0.0000296 |
4 |
0.0000039 |
0.5180040 |
| 512 |
NefscMQ2 |
s(PC1) |
0.8101715 |
4 |
1.0669234 |
0.0280692 |
| 512 |
NefscMQ2 |
s(PC2) |
0.8951363 |
4 |
2.1314256 |
0.0040081 |
| 512 |
NefscMQ2 |
s(PC3) |
0.1950601 |
4 |
0.0605743 |
0.2730047 |
| 513 |
NefscMQ2 |
s(PC1) |
0.8765322 |
4 |
1.7747800 |
0.0074476 |
| 513 |
NefscMQ2 |
s(PC2) |
0.9330253 |
4 |
3.4797996 |
0.0004663 |
| 513 |
NefscMQ2 |
s(PC3) |
0.0000028 |
4 |
0.0000002 |
0.6414642 |
Cluster Analysis
GOMindex_ca <- MEDMR_cpue %>%
filter(Season == "Spring") %>%
select(Year, lob_index, name, stat_area)
NEFSCindex_ca <- NEFSCindex %>%
filter(Season == 2) %>%
select(Year, lob_index, name)
lobData <- bind_rows(ALSI, sublegal_cpue, ME_pounds, GOMindex_ca, NEFSCindex_ca) %>%
pivot_wider(names_from = c(name, stat_area),
values_from = lob_index) %>%
na.omit() %>%
column_to_rownames("Year")
# standardize variables
lobData <- scale(lobData) %>%
data.frame()
lobData_t <- data.table::transpose(lobData)
# get row and colnames in order
colnames(lobData_t) <- rownames(lobData)
rownames(lobData_t) <- colnames(lobData)
# determine number of clusters
wss <- (nrow(lobData)-1)*sum(apply(lobData,2,var)) # get sum of squares
for (i in 2:12) wss[i] <- sum(kmeans(lobData,
centers=i)$withinss)
# look for break in plot like a scree plot
breakplot <- plot(1:12, wss, type="b", xlab="Number of Clusters",
ylab="Within groups sum of squares")

breakinfo <- fpc::pamk(lobData)
# K-Means Cluster Analysis
fit <- kmeans(lobData, 2) # 3 cluster solution
# get cluster means
cluster_means <- aggregate(lobData,by=list(fit$cluster),FUN=mean)
# append cluster assignment
lobData_cluster <- data.frame(lobData, fit$cluster)
# Cluster Plot against 1st 2 principal components
# vary parameters for most readable graph
cluster::clusplot(lobData, fit$cluster, color=TRUE, shade=TRUE,
labels=2, lines=0)

# Centroid Plot against 1st 2 discriminant functions
centroidPlot <- fpc::plotcluster(lobData, fit$cluster)
